#include "PO_MoM.h"
#include "RadiationPattern.h"
#include "FillMatrix.h"
#include <cmath>
#include <fstream>
using namespace blitz;
void Output()
{
	Pattern.resize(181);
	double max = 0;
	for (int theta=0; theta<=180; theta++)
	{
		Pattern(theta) = Farfld(theta/180.0*PI,0);
		if (Pattern(theta)>max)
			max = Pattern(theta);
	}
// 	Pattern = Pattern/max;
	std::ofstream fout("radiation_pattern.dat");
	for (int i=0; i<=180; i++)
		fout<<Pattern(i)<<std::endl;
	fout.close();
}
double Farfld(double theta, double phi)
{
	double thesin = sin(theta);
	double thecos = cos(theta);
	double phisin = sin(phi);
	double phicos = cos(phi);
	Position3D vRk(WaveNumber*thesin*phicos,WaveNumber*thesin*phisin,WaveNumber*thecos);
	Complex3D PO_Far(0,0,0);
	Complex3D MoM_Far(0,0,0);
	Complex3D Junction_Far(0,0,0);
	Complex3D Total_Far(0,0,0);
	for (int i=0; i<iEdge; i++)
	{
		int node1 = Edge(i,0);
		int node2 = Edge(i,1);
		Complex3D tmp(0.0,0.0,0.0);
		for (int LL=0; LL<2; LL++)
		{
			int pat = Edge(i,LL+2);
			int node3 = Patch(pat,0)+Patch(pat,1)+Patch(pat,2)-node1-node2;
			double signal(1.0-2*LL);

			Position3D vRho(vPatchCtr(pat)-vNode(node3));
			tmp +=vRho * signal * exp(dComplex(0.0,1.0)*dot(vRk,vPatchCtr(pat)));
		}
		PO_Far +=  Edge_length(i) * J(i) * tmp;
	}
	PO_Far = PO_Far * 0.5;
	for (int i=0; i<iWire; i++)
	{
 		MoM_Far += exp(dComplex(0.0,1.0) * dot(vRk,vWire(i))) * I(i);
	}
	MoM_Far = MoM_Far * Position3D(0.0,0.0,DeltaLength);
	for (int i=0; i<iCpatch; i++)
	{
		Junction_Far += exp(dComplex(0.0,1.0) * dot(vRk,vPatchCtr(Cpatch(i)))) * Areas(0,i) * f_k(i);
	}
	Junction_Far = Junction_Far * I(iWire-1);
	Total_Far = -dComplex(0.0,1.0)*Omega*Mu/(4*PI) * (MoM_Far+PO_Far+Junction_Far);
	return 20*log(abs(sqrt(dot(Total_Far,Total_Far))));
}
